home *** CD-ROM | disk | FTP | other *** search
/ Tech Arsenal 1 / Tech Arsenal (Arsenal Computer).ISO / tek-20 / ffteme12.zip / FFTEME12.PAS < prev    next >
Pascal/Delphi Source File  |  1991-04-02  |  11KB  |  498 lines
  1. {$R+}    {Range checking}
  2. {$B+}    {Boolean complete evaluation}
  3. {$S+}    {Stack checking}
  4. {$I+}    {I/O checking}
  5. {$N+}    {numeric coprocessor}
  6.  
  7. program FFTEME;
  8. { by Mike Cook (AF9Y) Mar 31, 1991}
  9.  
  10. Uses
  11.   DOS,
  12.   Crt,
  13.   Graph;
  14.  
  15. var
  16.   GraphDriver: integer;
  17.   GraphMode: integer;
  18.   ErrorCode: integer;
  19.  
  20. FileVar: File of Byte;
  21. FastFile: File;
  22. OutFile: Text;
  23. FileName: String[14];
  24.  
  25.   SampleByte: Byte;
  26.   ImagF: Array [0..4096] of Single;
  27.   RealF: Array [0..4096] of Single;
  28.   Sample: Array [0..28000] of Byte;
  29.  
  30.   Mag: Single;
  31.   Period: Single;
  32.  
  33.   TotalSamples: Integer;
  34.   TotalFileSamples: Integer;
  35.  
  36.   Ch: Char;
  37.   N: Integer;
  38.   NLimit: Integer;
  39.  
  40.   X1: Integer; X2: Integer;
  41.   Y1: Integer; Y2: Integer;
  42.  
  43.   I,X,L,Y,K,M,MX: Integer;
  44.  
  45.   LoopEnd: Integer;
  46.  
  47.   X3,Max: Single;
  48.   Z,Zmax: Integer;
  49.   ADwait: Integer;
  50.   GraphLine: Integer;
  51.   PreviousGraphLine: Integer;
  52.   YScale: Single;
  53.   YStart: Integer;
  54.   XOffSet: Integer;
  55.   YOffSet: Integer;
  56.  
  57.   NormFFTAns: Char;
  58.   NormFFT: Integer;
  59.   Start: Integer;
  60.   StringZ: String[40];
  61.   DefaultMinFreq, MinFreq: Integer;
  62.   DefaultMaxFreq, MaxFreq: Integer;
  63.   MinZ:    Integer;
  64.   MaxZ:    Integer;
  65.   DataSkips: Integer;
  66.   Skips: Integer;
  67.  
  68.   DisplayRow: Integer;
  69.   DisplayCol: Integer;
  70.  
  71.   HourA,MinuteA,SecondA,Sec100A: Word;
  72.   HourB,MinuteB,SecondB,Sec100B: Word;
  73.  
  74.   StartTime, EndTime: LongInt;
  75.  
  76.   Break: Boolean;
  77.  
  78. Procedure GraphSetup;
  79. Begin
  80.   GraphDriver:=Detect;
  81.   InitGraph(GraphDriver, GraphMode, '');
  82.   end;
  83.  
  84. Procedure InputSampleTotal;
  85. Var
  86.    OffsetLimit: Integer;
  87.  
  88. Begin
  89.   Write ('Number of Samples (4096,2048,1024)? '); Readln(TotalSamples);
  90.      L:=Round(Ln(TotalSamples)/Ln(2));
  91.      Period:=((TotalSamples/4096)*0.2485);{Recorded Sample Rate Cal aginist wwv}
  92.      ADwait:=Round(6+(28500*Period/TotalSamples));
  93.  
  94.   Write ('DataSkips (1,2,3,4)? '); Readln(DataSkips);
  95.      Period:=Period*DataSkips;
  96.      OffsetLimit:=Trunc(27776-(DataSkips*TotalSamples));
  97.  
  98.   Write ('Offset Samples (0..',OffsetLimit,')? '); Readln(Start);
  99.  
  100.   Write  ('Normal or Squared FFT (N/S)? '); Readln(NormFFTAns);
  101.      If ((NormFFTAns='S') or (NormFFTAns='s')) then NormFFT:=0 else NormFFT:=1;
  102.      If NormFFT=0 then Period:=Period*2;
  103.      DefaultMinFreq:=Round(1/Period);
  104.      DefaultMaxFreq:=Round(TotalSamples/2/Period);
  105.  
  106.   WriteLn('FFT Range = ',DefaultMinFreq,' Hz',' to ',DefaultMaxFreq,' Hz');
  107.  
  108.   Write ('Min FFT Freq? '); Readln(MinFreq);
  109.          MinZ:=Round(MinFreq*Period);
  110.  
  111.   Write ('Max FFT Freq? '); Readln(MaxFreq);
  112.          MaxZ:=Round(MaxFreq*Period);
  113.  
  114. End;
  115.  
  116. Procedure RequestReadFile;
  117. begin
  118.    Write('Enter name of file to read: ');
  119.    Readln(FileName);
  120.    {FileName:='EME1';}
  121.    Assign(FileVar,FileName);
  122.    Reset(FileVar);
  123. end;
  124.  
  125. Procedure FastReadFile;
  126. Var
  127.    NumRead: Word;
  128.  
  129. Begin
  130.    Assign(FastFile,FileName);
  131.    Reset(FastFile,1);
  132.  
  133.    Repeat
  134.       BlockRead(FastFile,Sample,SizeOf(Sample),NumRead);
  135.       Until (NumRead = 0);
  136.       Close(FileVar);
  137.  
  138.    NLimit:=TotalSamples-1;
  139.  
  140.  
  141. end;
  142.  
  143. Procedure ReadFile;
  144. Var
  145.    ReadCnt: Integer;
  146.    TotalReads: Single;
  147.    StoreCnt: Integer;
  148.  
  149. Begin
  150.      WriteLn;
  151.      DisplayRow:=WhereY;
  152.      DisplayCol:=WhereX;
  153.  
  154.      StoreCnt:=0; ReadCnt:=0; TotalReads:=(TotalSamples*DataSkips+Start)/100;
  155.  
  156.      While ((Not EOF(FileVar)) AND (StoreCnt<TotalSamples)) do
  157.         Begin
  158.            If ReadCnt<Start then
  159.               Begin
  160.                    Read(FileVar,SampleByte);
  161.                    Inc(ReadCnt);
  162.                    ClrEol; GotoXY(DisplayCol,DisplayRow);
  163.                    Write('Data Read ',(Round(ReadCnt/TotalReads)),'% Complete');
  164.                    end
  165.               else begin
  166.                    For Skips:=1 to DataSkips do
  167.                        Begin
  168.                        Read(FileVar,SampleByte);
  169.                        Inc(ReadCnt);
  170.                        If EOF(FileVar) then Skips:=DataSkips;
  171.                        end;
  172.                    Sample[(StoreCnt)]:=(SampleByte);
  173.                    StoreCnt:=StoreCnt+1;
  174.                    ClrEol; GotoXY(DisplayCol,DisplayRow);
  175.                    Write('Data Read ',(Round(ReadCnt/TotalReads)),'% Complete');
  176.                    end;
  177.          end;
  178.      Close(FileVar);
  179.      If StoreCnt<TotalSamples then
  180.         Begin
  181.             Writeln('Not Enough Data');
  182.             Halt;
  183.             end;
  184. End;
  185.  
  186. Procedure DisableClock;
  187. Begin
  188.    Port[$21]:=(Port[$21] or $01);
  189. End;
  190.  
  191. Procedure EnableClock;
  192. Begin
  193.    Port[$21]:=(Port[$21] and $FE);
  194. End;
  195.  
  196. Procedure GenSineWave;
  197. Var
  198.   Freq: Single;
  199.   TimeStep: Single;
  200.   TimeN: Single;
  201.   Period: Single;
  202.  
  203. Begin
  204.   Freq:=555; TimeStep:=Period/Totalsamples;
  205.   TimeN:=0;
  206.   NLimit:=TotalSamples-1;
  207.   For N:=0 to Nlimit do
  208.   Begin
  209.     Sample[N]:=Round(128*(Sin(2*PI*Freq*TimeN)));
  210.     TimeN:=TimeN+TimeStep;
  211.   End;
  212. End;
  213.  
  214. Procedure CopyFastSamples;
  215. Var
  216.   IndexN: Integer;
  217.  
  218. Begin
  219.   NLimit:=TotalSamples-1;
  220.   IndexN:=Start;
  221.   For N:=0 to NLimit do
  222.   Begin
  223.      If NormFFT=1 then RealF[N]:=((Sample[IndexN]-128)/128)/TotalSamples
  224.                        { Normal FFT }
  225.         else begin
  226.         RealF[N]:=(((Sample[IndexN]-128)*(Sample[IndexN]-128))/128-128)/TotalSamples;
  227.                        { Square FFT }
  228.         end;
  229.  
  230.      ImagF[N]:=0;
  231.      IndexN:=IndexN+DataSkips;
  232.      end;
  233. End;
  234.  
  235. Procedure CopySamples;
  236. Begin
  237.   NLimit:=TotalSamples-1;
  238.   For N:=0 to NLimit do
  239.   Begin
  240.      If NormFFT=1 then RealF[N]:=(Sample[N]/128-128)/TotalSamples      { Normal FFT }
  241.         else
  242.      RealF[N]:=((Sample[N])*Sample[N]/128-128)/TotalSamples;         { Square FFT }
  243.      ImagF[N]:=0;
  244.   end;
  245. End;
  246.  
  247. Procedure Scramble;
  248. Var
  249.   W: Integer;
  250.   N1: Integer;
  251.  
  252. Begin
  253.   Y:=0; N1:=TotalSamples;
  254.   For W:=1 to L do
  255.   Begin
  256.     N1:=N1 shr 1;
  257.     If X>=N1 then
  258.     Begin
  259.       If W=1 then Inc(Y)
  260.       else
  261.       Y:=Y+(2 shl (W-2));
  262.       X:=X-N1;
  263.     end;
  264.   end;
  265. End;
  266.  
  267. Procedure CalFFT;
  268. Var
  269.    I1,I2,I3,I4,I5: Integer;
  270.    V,A1,A2,B1,B2,Z1,Z2: Single;
  271.  
  272. Begin
  273. I1:=TotalSamples shr 1; I2:=1; V:=2*Pi/TotalSamples;
  274. WriteLn;
  275.  DisplayRow:=WhereY;
  276.  DisplayCol:=WhereX;
  277.  
  278.   For I:=1 to L do
  279.   Begin
  280.     ClrEol; GotoXY(DisplayCol,DisplayRow);
  281.     Write('FFT ',(Round((I/L)*100)),'% Complete');
  282.  
  283.     I3:=0; I4:=I1;
  284.  
  285.     For K:=1 to I2 do
  286.     Begin
  287.  
  288.        X:=Trunc(I3/I1); Scramble;
  289.        I5:=Y;
  290.        Z1:=Cos(V*I5); Z2:=-Sin(V*I5);
  291.  
  292.        LoopEnd:=I4-1;
  293.        For M:=I3 to LoopEnd do
  294.        Begin
  295.           A1:=RealF[M]; A2:=ImagF[M]; Mx:=M+I1;
  296.           B1:=Z1*RealF[Mx]-Z2*ImagF[Mx];
  297.           B2:=Z2*RealF[Mx]+Z1*ImagF[Mx];
  298.           RealF[M]:=(A1+B1); ImagF[M]:=(A2+B2);
  299.           RealF[Mx]:=(A1-B1); ImagF[Mx]:=(A2-B2);
  300.        end;
  301.  
  302.        I3:=I3+(I1 shl 1); I4:=I4+(I1 shl 1);
  303.     end;
  304.     I1:=I1 shr 1; I2:=I2 shl 1;
  305.   end;
  306. end;
  307.  
  308. Procedure CalMagnitude;
  309. Begin
  310.      Scramble;
  311.      Mag:=SQRT(RealF[Y]*RealF[Y]+ImagF[Y]*ImagF[Y]);
  312.  
  313. End;
  314.  
  315. Procedure CalMaxSpectrum;
  316. Begin
  317.      Max:=0;
  318.      LoopEnd:=TotalSamples shr 1;
  319.      For Z:=MinZ to MaxZ do
  320.      Begin
  321.         X:=Z;
  322.         CalMagnitude;
  323.         If Mag>Max then
  324.          Begin
  325.             Max:=Mag;
  326.             Zmax:=Z;
  327.             End;
  328.      End;
  329. End;
  330.  
  331. Procedure PlotTimeSamples;
  332. Begin
  333.      YScale:=0.13;  YStart:=10;
  334.      GraphLine:=0; PreviousGraphLine:=0;
  335.      LoopEnd:=TotalSamples-1;
  336.  
  337.      For N:=0 to LoopEnd do
  338.      Begin
  339.         Case N of
  340.             0000..0639 : GraphLine:=1;
  341.             0640..1279 : GraphLine:=2;
  342.             1280..1919 : GraphLine:=3;
  343.             1920..2559 : GraphLine:=4;
  344.             2560..3199 : GraphLine:=5;
  345.             3200..3839 : GraphLine:=6;
  346.             3840..4479 : GraphLine:=7;
  347.             4480..5120 : GraphLine:=8;
  348.             end;
  349.  
  350.         If Graphline<>PreviousGraphLine then
  351.             Begin
  352.                YOffSet:=YStart+GraphLine*35;
  353.                Y1:=Round(YOffSet-YScale*(Sample[N]));
  354.                XOffSet:=Round((GraphLine-1)*640);
  355.                PreviousGraphLine:=GraphLine;
  356.                end;
  357.  
  358.         Y2:=Round(YOffSet-YScale*(Sample[N]+128));
  359.         Line(N-XOffSet,Y1,N-XOffSet,Y2);
  360.         Y1:=Y2
  361.         End;
  362. end;
  363.  
  364. Procedure PlotFastTimeSamples;
  365. Var
  366.      IndexN: Integer;
  367.  
  368. Begin
  369.      YScale:=0.13;  YStart:=10;
  370.      GraphLine:=0; PreviousGraphLine:=0;
  371.      LoopEnd:=TotalSamples-1;
  372.  
  373.      For N:=0 to LoopEnd do
  374.      Begin
  375.         Case N of
  376.             0000..0639 : GraphLine:=1;
  377.             0640..1279 : GraphLine:=2;
  378.             1280..1919 : GraphLine:=3;
  379.             1920..2559 : GraphLine:=4;
  380.             2560..3199 : GraphLine:=5;
  381.             3200..3839 : GraphLine:=6;
  382.             3840..4479 : GraphLine:=7;
  383.             4480..5120 : GraphLine:=8;
  384.             end;
  385.  
  386.         IndexN:=N*DataSkips+Start;
  387.         If Graphline<>PreviousGraphLine then
  388.             Begin
  389.                YOffSet:=YStart+GraphLine*35;
  390.                Y1:=Round(YOffSet-YScale*(Sample[IndexN]));
  391.                XOffSet:=Round((GraphLine-1)*640);
  392.                PreviousGraphLine:=GraphLine;
  393.                end;
  394.  
  395.         Y2:=Round(YOffSet-YScale*(Sample[IndexN]+128));
  396.         Line(N-XOffSet,Y1,N-XOffSet,Y2);
  397.         Y1:=Y2
  398.         End;
  399. end;
  400.  
  401. Procedure PlotSpectrum;
  402. Begin
  403.      YScale:=0.40;  YStart:=260;
  404.      GraphLine:=0; PreviousGraphLine:=0;
  405.      LoopEnd:=TotalSamples Shr 1;
  406.  
  407.      For Z:=MinZ to MaxZ do
  408.      Begin
  409.         Case Z of
  410.             0000..0639 : GraphLine:=1;
  411.             0640..1279 : GraphLine:=2;
  412.             1280..1919 : GraphLine:=3;
  413.             1920..2559 : GraphLine:=4;
  414.             2560..3200 : GraphLine:=5;
  415.             3200..3839 : GraphLine:=6;
  416.             3840..4480 : GraphLine:=7;
  417.             4480..5120 : GraphLine:=8;
  418.             end;
  419.  
  420.         If Graphline<>PreviousGraphLine then
  421.             Begin
  422.                YOffSet:=Round(YStart+GraphLine*50.0);
  423.                XOffSet:=Round((GraphLine-1)*640.0);
  424.                PreviousGraphLine:=GraphLine;
  425.                end;
  426.         X:=Z;
  427.         CalMagnitude;
  428.         Y1:=Round(YOffSet-(YScale*Mag/Max)*199.0);
  429.         Line(Z-XOffSet,YOffSet,Z-XOffSet,Y1);
  430.         End;
  431. end;
  432.  
  433. Procedure PrintSpectrum;
  434. Begin
  435.      writeln;
  436.      LoopEnd:=TotalSamples shr 1;
  437.      For Z:=0 to LoopEnd do
  438.      Begin
  439.         X:=Z;
  440.         CalMagnitude;
  441.         WriteLn(Z,' ',RealF[Y],' ',ImagF[Y],' ',Mag);
  442.      End;
  443. end;
  444.  
  445. {=============================  MAIN ====================================}
  446.  
  447.  
  448. Begin
  449.  
  450. Break:=False;
  451.  
  452. Repeat
  453.  
  454. Begin
  455. TextMode(CO80);
  456.  
  457. ClrScr;
  458.  
  459. InputSampleTotal;
  460. {GenSineWave;}
  461.  
  462. RequestReadFile;
  463.    {ReadFile;}
  464.    FastReadFile;
  465.  
  466.    {CopySamples;}
  467.    CopyFastSamples;
  468.  
  469. GetTime(HourA,MinuteA,SecondA,Sec100A);
  470. CalFFT;
  471.        GetTime(HourB,MinuteB,SecondB,Sec100B);
  472.        StartTime:=HourA*360000+MinuteA*6000+SecondA*100+Sec100A;
  473.        EndTime:=HourB*360000+MinuteB*6000+SecondB*100+Sec100B;
  474.        Writeln;
  475.        Writeln('FFT Process Time = ',((EndTime-StartTime)/100):2:2);
  476.  
  477. CalMaxSpectrum;
  478.  
  479. GraphSetup;
  480.  
  481. PlotFastTimeSamples;
  482. {PlotTimeSamples;}
  483.  
  484. {PrintSpectrum;}
  485. PlotSpectrum;
  486.  
  487. Zmax:=Round(Zmax/Period);
  488.       Str(Zmax,StringZ); OutTextXY(600,460,StringZ);
  489.  
  490. Repeat Until Keypressed;
  491.  
  492. end;
  493.  
  494. Until Break;
  495.  
  496. TextMode(CO80);
  497. End.
  498.